!***********************************************************
!
program inc_to_tec
!BE CAREFUL WITH FORTRAN COMPILER
!USE OPTION ifort -assume byterecl
!***********************************************************

implicit none



! Variables reading from Incompact3D, real(8) for double precision
integer(4),parameter :: nx=513,ny=193,nz=256,xjicf=372,ibm_y_max=97
real(8), parameter :: lx=24., ly=4., lz=6.

include "declare.f90"
include "tecio.f90"

l=1
990 format('phi',I3.3)
write(nf,990) l
inquire(file=trim(nf), exist=iscalar)

if (.not. iscalar) then
  Title="JICF - 3D - no-recycling - no Phi - explicit - instantaneous"
  Variables="X Y Z p u v w"
  n_variable = 4
else
  Title="JICF - 3D - no-recycling - high sc - explicit - instantaneous"
  Variables="X Y Z p u v w phi"
  n_variable = 5
endif
FName="ins_JICF.plt"

!write(*,*) n_variable,' files to read and export for Tecplot'
allocate (data(nx,ny,nz,n_variable+3))
allocate (epsilon(nx,ny,nz))
allocate (PassiveVarList(n_variable+3))
allocate (ValueLocation(n_variable+3))
allocate (ShareVarFromZone(n_variable+3))
PassiveVarList=0
ValueLocation=1 !ValueLocation: 0=cell-centered, 1=node-centered
ShareVarFromZone=0
IJKN = nx*ny*nz*(n_variable+3)

!dx and dz
!dx=lx/nx !periodic
dx=lx/(nx-1) !non-periodic
dz=lz/nz !periodic
!dz=lz/(nz-1) !non-periodic
 
! Reading the mesh
write(*,*) 'Reading the grids y'
open(21,file='yp.dat',form='formatted')
do j=1,ny
  read(21,*) i,yp(j)
enddo
close(21)
write(*,*) 'Grids y are read'

! reading IBM
if (.false.) then
  OPEN(21,FILE=trim('epsilon.dat'),FORM='UNFORMATTED',&
	ACCESS='DIRECT', RECL=8, STATUS='OLD')!, CONVERT='big_endian')
  data_count = 1
  DO k=1,nz
  DO j=1,ny
  DO i=1,nx
    READ(21,REC=data_count) epsilon(i,j,k)
    data_count = data_count + 1
  ENDDO
  ENDDO
  ENDDO
  CLOSE(21)
else
  DO k=1,nz
  DO j=1,ny
  DO i=1,nx
    tempa=xjicf-i
    tempb=nz/2-k
    r2=dx*dx*real(tempa*tempa)+dz*dz*real(tempb*tempb)
    if ((j .le. ibm_y_max) .and. (r2 > 1.)) then
      epsilon(i,j,k)=1.
    else
      epsilon(i,j,k)=0.
    endif
  ENDDO
  ENDDO
  ENDDO
endif
write(*,*) 'Reading IBM completed'
!epsilon=0.


! Calculating coordinates for Paraview
DO k=1,nz
DO j=1,ny
DO i=1,nx
  data(i,j,k,1)=real(   (0+(i-xjicf))*dx) !be careful about (0+(i-1)*step) or (1+(i-1)*step)
  data(i,j,k,2)=real( yp(1+(j-1))   )
  data(i,j,k,3)=real(   (0+(k-nz/2))*dz) !be careful about (0+(k-1)*step) or (1+(k-1)*step)
ENDDO
ENDDO
ENDDO

err = tecini112 &
  (trim(Title)//char(0), &
  trim(Variables)//char(0), &
  trim(FName)//char(0), &
  "."//char(0), &
  FileType, &
  Debug, &
  VIsDouble)
write(*,*) '-----------Initializing completed-----------' 

123 format('pp',I3.3)
write(nf,123) l
write(*,*) 'Finding file ',trim(nf)
inquire(file=trim(nf), exist=dir_e)

do while (dir_e)

!****************************************************
! Pressure
!****************************************************
  ! reading file
991 format('pp',I3.3)
  write(nf,991) l
  !write(*,*) 'Reading file ',trim(nf)
  ! ATTENTION: SI DATA EN PROVENANCE IDRIS ALORS CONVERT='big_endian' (IF DATA FROM IDRIS THEN CONVERT)
  OPEN(21,FILE=trim(nf),FORM='UNFORMATTED',&
    ACCESS='DIRECT', RECL=8, STATUS='OLD')!, CONVERT='big_endian')
  data_count = 1
  DO k=1,nz
  DO j=1,ny
  DO i=1,nx
    READ(21,REC=data_count) data(i,j,k,4)
    if (epsilon(i,j,k)==1.) data(i,j,k,4)=0.
    data_count = data_count + 1
  ENDDO
  ENDDO
  ENDDO
  CLOSE(21)
  write(*,*) 'Reading file ',trim(nf),' completed'
 
!****************************************************
! Ux
!****************************************************
  ! reading file
992 format('ux',I3.3)
  write(nf,992) l
  !write(*,*) 'Reading file ',trim(nf)
  OPEN(21,FILE=trim(nf),FORM='UNFORMATTED',&
    ACCESS='DIRECT', RECL=8, STATUS='OLD')!, CONVERT='big_endian')
  data_count = 1
  DO k=1,nz
  DO j=1,ny
  DO i=1,nx
    READ(21,REC=data_count) data(i,j,k,5)
    if (epsilon(i,j,k)==1.) data(i,j,k,5)=0.
    data_count = data_count + 1
  ENDDO
  ENDDO
  ENDDO
  CLOSE(21)
  write(*,*) 'Reading file ',trim(nf),' completed'
 
!****************************************************
! Uy
!****************************************************
  ! reading file
993 format('uy',I3.3)
  write(nf,993) l
  !write(*,*) 'Reading file ',trim(nf)
  OPEN(21,FILE=trim(nf),FORM='UNFORMATTED',&
    ACCESS='DIRECT', RECL=8, STATUS='OLD')!, CONVERT='big_endian')
  data_count = 1
  DO k=1,nz
  DO j=1,ny
  DO i=1,nx
    READ(21,REC=data_count) data(i,j,k,6)
    if (epsilon(i,j,k)==1.) data(i,j,k,6)=0.
    data_count = data_count + 1
  ENDDO
  ENDDO
  ENDDO
  CLOSE(21)
  write(*,*) 'Reading file ',trim(nf),' completed'

!****************************************************
! Uz
!****************************************************
  ! reading file
994 format('uz',I3.3)
  write(nf,994) l
  !write(*,*) 'Reading file ',trim(nf)
  OPEN(21,FILE=trim(nf),FORM='UNFORMATTED',&
    ACCESS='DIRECT', RECL=8, STATUS='OLD')!, CONVERT='big_endian')
  data_count = 1
  DO k=1,nz
  DO j=1,ny
  DO i=1,nx
    READ(21,REC=data_count) data(i,j,k,7)
    if (epsilon(i,j,k)==1.) data(i,j,k,7)=0.
    data_count = data_count + 1
  ENDDO
  ENDDO
  ENDDO
  CLOSE(21)
  write(*,*) 'Reading file ',trim(nf),' completed' 

!****************************************************
! Phi
!****************************************************
  ! reading file
  if (iscalar) then
  write(nf,990) l
  !write(*,*) 'Reading file ',trim(nf)
  OPEN(21,FILE=trim(nf),FORM='UNFORMATTED',&
    ACCESS='DIRECT', RECL=8, STATUS='OLD')!, CONVERT='big_endian')
  data_count = 1
  DO k=1,nz
  DO j=1,ny
  DO i=1,nx
    READ(21,REC=data_count) data(i,j,k,8)
    if (epsilon(i,j,k)==1.) data(i,j,k,8)=0.
    data_count = data_count + 1
  ENDDO
  ENDDO
  ENDDO
  CLOSE(21)
  write(*,*) 'Reading file ',trim(nf),' completed'
  endif

  SolTime=l

  !/*  Ordered Zone */
  err = teczne112 &
    ("Ordered Zone"//char(0), &
    ZoneType, &
    nx, &
    ny, &
    nz, &
    ICellMax, &
    JCellMax, &
    KCellMax, &
    SolTime, &
    StrandID, &
    ParentZn, &
    IsBlock, &
    NFConns, &
    FNMode, &
    TotalNumFaceNodes, &
    TotalNumBndryFaces, &
    TotalNumBndryConnections, &
    PassiveVarList, &
    ValueLocation, &
    ShareVarFromZone, &
    ShrConn)

  err = tecdat112(IJKN, data, DIsDouble)

800 format('Time step ',I3.3)
  write(nf,800) l
  write(*,*) '-----------',trim(nf),' completed-----------' 

  l=l+1
  write(nf,123) l
  inquire(file=trim(nf), exist=dir_e)

enddo
! end loop

write(*,*) '-----------Can''t find ',trim(nf),'-----------'
err = tecend112();

write(*,*) '-----------END-----------'

end program inc_to_tec